//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2012 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : matv_kernels.h
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#ifndef MATV_KERNELS_H
#define	MATV_KERNELS_H

namespace rrbsolver
{

template <class T>
__global__ void kernel_matv1(T *y, const T *x, const Grid g)
{
    int ld = g.nx;

    int u_r1 = BORDER_WIDTH + bx * Bx + tx;
    int v_r1 = BORDER_WIDTH + by * By + ty;
    int u_r2 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_r2 = BORDER_WIDTH2 + g.cy + by * By + ty;

    int u_b1 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_b1 = BORDER_WIDTH + by * By + ty;
    int u_b2 = BORDER_WIDTH + bx * Bx + tx;
    int v_b2 = BORDER_WIDTH2 + g.cy + by * By + ty;

    int loc_r1 = ld * v_r1 + u_r1;
    int loc_r2 = ld * v_r2 + u_r2;
    int loc_b1 = ld * v_b1 + u_b1;
    int loc_b2 = ld * v_b2 + u_b2;

    const T x_r1    = __FETCH_XX(u_r1,   v_r1,   ld);
    const T x_r1up1 = __FETCH_XX(u_r1+1, v_r1,   ld);
    const T x_r1vp1 = __FETCH_XX(u_r1,   v_r1+1, ld);
    const T x_r2    = __FETCH_XX(u_r2,   v_r2,   ld);
    const T x_r2vm1 = __FETCH_XX(u_r2,   v_r2-1, ld);
    const T x_r2um1 = __FETCH_XX(u_r2-1, v_r2,   ld);

    const T nn_r1   = g.nn[loc_r1];
    const T ee_r1   = g.ee[loc_r1];
    const T ss_r2   = g.ss[loc_r2];
    const T ww_r2   = g.ww[loc_r2];

    // b1:
    T sum_b1 = g.cc[loc_b1] * (ss_r2        * x_r2 +     // north
                               g.ee[loc_b1] * x_r1up1 +  // east
                               g.ss[loc_b1] * x_r2vm1 +  // south
                               ee_r1        * x_r1);     // west

    // b2:
    T sum_b2 = g.cc[loc_b2] * (g.nn[loc_b2] * x_r1vp1 +  // north
                               ww_r2        * x_r2 +     // east
                               nn_r1        * x_r1 +     // south
                               g.ww[loc_b2] * x_r2um1);  // west

    y[loc_b1] = sum_b1;
    y[loc_b2] = sum_b2;

    // r1 (partial):
    T sum_r1 = g.cc[loc_r1] * x_r1 -
               nn_r1        * sum_b2 -
               ee_r1        * sum_b1;


    // r2 (partial):
    T sum_r2 = g.cc[loc_r2] * x_r2 -
               ss_r2        * sum_b1 -
               ww_r2        * sum_b2;

    y[loc_r1] = sum_r1;
    y[loc_r2] = sum_r2;
}

template <class T>
__global__ void kernel_matv2(T *y, const T *x, const Grid g)
{
    int ld = g.nx;

    int u_r1 = BORDER_WIDTH + bx * Bx + tx;
    int v_r1 = BORDER_WIDTH + by * By + ty;
    int u_r2 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_r2 = BORDER_WIDTH2 + g.cy + by * By + ty;

    int u_b1 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_b1 = BORDER_WIDTH + by * By + ty;
    int u_b2 = BORDER_WIDTH + bx * Bx + tx;
    int v_b2 = BORDER_WIDTH2 + g.cy + by * By + ty;

    int loc_r1 = ld * v_r1 + u_r1;
    int loc_r2 = ld * v_r2 + u_r2;

    const T y_b2vm1 = __FETCH_YY(u_b2,   v_b2-1, ld);
    const T y_b1um1 = __FETCH_YY(u_b1-1, v_b1,   ld);
    const T y_b1vp1 = __FETCH_YY(u_b1,   v_b1+1, ld);
    const T y_b2up1 = __FETCH_YY(u_b2+1, v_b2,   ld);

    // r1 (correction):
    T sum_r1 = y[loc_r1] - g.ss[loc_r1] * y_b2vm1  // south
                         - g.ww[loc_r1] * y_b1um1; // west

    // r2 (correction):
    T sum_r2 = y[loc_r2] - g.nn[loc_r2] * y_b1vp1  // north
                         - g.ee[loc_r2] * y_b2up1; // east

    y[loc_r1] = sum_r1;
    y[loc_r2] = sum_r2;
}

} // namespace rrbsolver

#endif // MATV_KERNELS_H

